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Abstract 



Tissot indicatrices have provided visual measures of local area and isotropy distortions. 
Here we show how large scale distortions of flexion (bending) and skewness (lopsidedness) can 
be measured. Area and isotropy distortions depend on the map projection metric, flexion and 
skewness, which manifest themselves on continental scales, depend on the first derivatives of 
the metric. We introduce new indicatrices that show not only area and isotropy distortions 
but flexion and skewness as well. We present a table showing error measures for area, isotropy, 
flexion, skewness, distances, and boundary cuts allowing us to compare different world map 
projections. We find that the Winkel-Tripel projection (already adopted for world maps by the 



National Geographic), has low distortion on most measures and excellent quality overall. 
Keywords: Maps, Earth, Projection, Curvature 
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£j . 1 Introduction 

CO 

Tissot (1881) indicatrices have been very useful for providing a visual presentation of local dis- 
tortions in map projections in a simple and compelling fashion. A small circle of tiny radius (say 
0.1 degree of arc in radius) is constructed at a given location, and then enlarged and projected 
on the map at that location. This always produces an ellipse. Usually one favors conformal map 
projections that minimize the changes in scale factor, or equal area projections that minimize 
anisotropy, or recently, map projections that are neither conformal nor equal area, but which have 
a judicious combination of minimizing both scale and isotropy errors (like the Winkel-tripel used 
by the National Geographic Society for world maps). 

The Tissot ellipse at a given location is specified by three parameters, the major axis, the minor 
axis, and the orientation angle 9 of the major axis of the ellipse. Geometrically, from differential 
geometry (and General Relativity) we know that the measurement of local distances is measured 
by the metric tensor g a b, where a and b can each take the values x or y, yielding three independent 
components. Locally for two nearby points separated by infinitesimal map coordinate differences 
dx and dy, the true distance between these two points on the globe is given by: 

ds 2 = g xx dx 2 + 2g xy dx dy + g yy dy 2 (1) 

The Tissot ellipse (major axis, minor axis, and orientation angle, 6) can be calculated from the 
three components of the metric tensor. Thus, the Tissot ellipse essentially carries the information 
on the metric tensor for the map. It tells us how local infinitesimal distances on the map correlate 
with local infinitesimal distances on the globe. 
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As we will see, however, the Tissot ellipse does not carry all of the information related to 
distortions, which has been noted previously. Others (Stewart 1943; Peters 1975, 1978; Albinus 
1981; Canters 1989, 2002; Laskowski 1997ab) have noted that there are finite distortions apart from 
those of the Tissot ellipse. 

Previous authors have previously depicted finite size distortions by several methods, includ- 
ing showing faces on a globe (Reeves 1910; Gedymin 1946), 30 degree x 30 degree equiangular 
quadralaterals (Chamberlin 1947), a net of twenty spherical triangles, in an icosahedral arrange- 
ment (Fisher & Miller 1944), 150 degree great circle arcs, and circles of 15 degrees radius Tobler 
(1964). 

In the Oxford Hammond Atlas of the World (1993) , new conformal map projections ("Ham- 
mond Optimal Conformal Projections") were designed for the continents. Following the Chebyshev 
criterion (Snyder 1993), the rms scale factor errors were minimized by producing a constant scale 
factor along the boundary of the continent. (For a circular region this conformal map would be 
a stereographic projection.) By tailoring the boundary to the shape of each continent, the errors 
could be reduced relative to those in a simple stereographic projection. In touting the advantages 
of their projection the Hammond Atlas did the following experiment. Following Reeves (Reeves 
1910), they constructed a face on the globe with a triangular nose, a straight (geodesic) mouth, 
and eyes that were pairs of concentric circles on the globe. They then showed this face with various 
map projections. In the gnomonic, the mouth was straight, but the eye circles were not circular 
and were not concentric (what we would call skewness). The Mercator projection had the mouth 
smiling (what we would call flexion) and although the eyes were circular they were not concentric. 
The Hammond Optimal Conformal projections did a bit better on these qualities because the gra- 
dients of the scale factor changes were small so the flexion and skewness were small, although of 
course not zero. 

In Sec. [2] we introduce the concept of "flexion", by which a map projection can cause artificial 
bending of large structures. In Sec. [31 we show another distortion: skewness, which represents 
lopsidedness and an asymmetric stretching of large structures. We show a simple way to visualize 
these distortions in Sec. HI in which we introduce the "Goldberg-Gott Indicatrices." In Sec. [5] we 
derive a differential geometry approach to measuring flexion and skewness. While readers interested 
in computing flexion and skewness on projections not included in this paper should refer to Sec. 
it is highly technical, and those interested in seeing results may skip directly to Sec. El in which 
we discuss Monte Carlo estimates of the distortions for a number of projections, a ranking of map 
projections, and our conclusions. 



2 First Distortion: Flexion 

The local effects shown by the Tissot Ellipse are not the only distortions present in maps. There 
are also "flexion" (or bending) and "skewness" (or lopsidedness; discussed in the next section), 
which describe curvature distortions visible on world maps (the terminology stems from a similar 
effect in gravitational lensing; Goldberg & Bacon 2005). 

One can think about flexion in the following way. Imagine a truck going along a geodesic of 
the globe at unit angular speed (say one radian per day). Now imagine the image of that truck 
on the map, moving along. If the map were perfect, if it had zero flexion and zero skewness, then 
that truck would move in a straight line on the map with constant speed. Its velocity vector on 
the map: 

dx , , 
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would be a constant, where r is the angle of arclength in radians traveled by the truck along the 
geodesic on the globe. Its acceleration: 

would be zero. Of course, this cannot be true for a general geodesic. In the general case, the 
image of the truck suffers an acceleration as it moves along. The acceleration vector, a, in the 
two-dimensional map has two independent components: aj_ (which is perpendicular to the truck's 
velocity vector at that point), and an (which is parallel to its velocity vector at that point). 

The perpendicular acceleration, a±, causes the truck to turn without changing its speed on the 
map. This causes flexion, or bending, of geodesies. We define the flexion along a given geodesic at 
a given point to be: 

or, more usefully: 

/=^i. (5) 

GST V 

In this form, we may define: 

da = — — . (6) 

v 

where a is the angle of rotation suffered by the velocity vector. Remember, if a± is the only 
acceleration present, then the velocity vector of the truck on the map does not change in magnitude, 
but just rotates in angle. Thus / = da/dr, and represents the angular rate at which the velocity 
vector rotates divided by the angular rate at which the truck moves on the globe. 

Skewness and flexion only express themselves on large scales. They are not noticeable on 
infinitesimal scales where the metric contains all the information one needs, but become noticeable 
on finite scales, with their importance growing with the size of the area being examined. Flexion 
and skewness are thus important on continental scales and larger in a world map. 

It is possible to design a map projection that has zero flexion: the gnomonic projection shows 
all great circles as straight lines. However it does exhibit anisotropy, scale changes, and skewness, 
and at best can show only one hemisphere of the globe. 



2.1 Example 1: The Stereographic Projection 

As an example, consider a truck traveling on the equator as seen by in the polar stereographic 
projection (see Fig[T]). In the stereographic projection, the north pole is in the center of the map 
and the equator is a circle around it. As the truck circles the equator (the equator is a geodesic-so 
the truck drives straight ahead on the globe), it travels around a circle on the map. By azimuthal 
symmetry, the truck circles the equatorial circle on the map at a uniform rate. The velocity vector 
of the truck on the map rotates a complete 360° (2tt radians), as the truck circles the equator, 
traversing 360° of arc on the globe. So the flexion is f = 1, for a point on the equator, for a geodesic 
pointing in the direction of the equator. The flexion is defined at a point, and for a specific geodesic 
traveling through that point. 

For an arbitrary point on the equator in the stereographic projection and a geodesic pointing in 
the direction of a meridian of longitude (also a geodesic) the flexion is zero, because these geodesies 
are shown as radial straight lines in the polar stereographic projection, and the velocity vector of 
the truck does not turn as it travels north. 

The stereographic projection has the property that every great circle (geodesic) on the globe is 
shown as a circle on the map, except for a set of measure zero that pass through the north pole 
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(i.e. the meridians of longitude). Thus by the argument given above, the average flexion integrated 
around a random great circle must be (/) = 1, because the truck's velocity vector on a random 
great circle must rotate by 360° as it circles the 360° of arc completing that great circle on the 
globe. The magnitude of the velocity vector of the truck on the map is larger the further from the 
north pole it is and so its rotation per angle of arc of truck travel on the globe is larger there as 
well, and so the flexion along that random geodesic is larger the further away from the pole one is, 
with the integrated average value along the whole great circle being (/) = 1. 

2.2 Example 2: The Mercator Projection 

The Mercator projection (see Fig [2]) is conformal and so only the scale factor changes as a function 
of position on the map (i.e. g xx = g yy , and g xy = and the Tissot ellipses are all circles with radii 
proportional to l/g xx )- But there is bending. The northern boundary between the continental 
United States and Canada at the 49 th parallel of latitude is shown as a straight line in the Mercator 
Map, but really it is a small circle that is concave to the north. If one drove a truck down that 
border from west to east, one would have to turn the steering wheel slightly to the left so that one 
was continually changing direction. The great circle route (the straightest route) connecting the 
Washington State and Minnesota (both at the 49th parallel) is a straight line which goes entirely 
through Canada. This straight line on the globe when extended, passes south of the northern part 
of Maine, so the continental United States is bend downward like a frown in the Mercator Map. 
(See Figure [3] and H]). Likewise, Maine on a Mercator map Maine sags below the line connecting 
Washington State and Minnesota, while on the globe this is not true. 

In the Mercator projection, the flexion along the equator is zero, also along all meridians of 
longitude, but these are a set of measure zero. A random geodesic is a great circle that is inclined 
at some angle between 0° and 90° with respect to the equator. On the Mercator map this is a wavy 
line that bends downward in the northern hemisphere, and by symmetry, upward in the southern. 
Since the curvatures are equal and opposite in the two hemispheres, the average flexion (/) = 0, 
but this is misleading because the flexion at each point off the equator is not zero. So if we are 
rating map projections by the amount of flexion they contain we should use the absolute value of 
the flexion instead: |/|. In a region where the flexion does not change sign (such as the northern 
hemisphere in the Mercator projection or the entire stereographic map) the total bending of a 
geodesic segment will be the integral of the flexion |/| over that segment. In fact, in Section [6.11 
we will evaluate the overall flexion on a map by simply picking random points on the sphere and 
random directions for geodesies going through them, and then calculating the absolute value for 
the flexion for all random points on the globe and random directions through them. 

2.3 A Global Flexion Measure 

We can calculate the flexion for any point in the Mercator (or any other) projection through any 
geodesic using spherical trigonometry. As a reminder to the reader, the Mercator projection uses 
the mapping: 

x = A (7) 
y = ln(tan[7r/4 + 0/2]) 

(8) 

where here and throughout, A is the longitude expressed in radians, and (ft is the latitude expressed 
in radians. 
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On the map, the angle rotation along the geodesic with an azimuth, 9, is calculated by con- 
structing thin spherical triangle with side-angle-side given by (it/ 2 — 4>,6,dr). 

da = 7T-f3-9 (9) 

because the geodesic intersects the north-south meridian (a vertical line in the Mercator) at an 
angle of 9 initially, and at an angle tt — (3 at the other end, where (3 is the angle in the spherical 
triangle at the other end of the dr side. Solving for da using spherical trigonometry in the limit as 
dr goes to zero, we find that 

/ = sinOtancj) (10) 
Thus, for 9 = tt/2, an east-west geodesic, we find that 

f E w = tan<f> (11) 

so that in the northern hemisphere, traveling east one's geodesic is bending clockwise with da/dr = 
tan<f>. Therefore, the east-west geodesic bends downward. For, 9 = 0, a north-south geodesic, the 
flexion is zero, as we expect, since the meridians of longitude are straight in the Mercator map. If 
we average over all azimuths at a given point, we find: 

<l/(0)l> = |tan0|- (12) 

TT 

Now we can integrate this over all points on the sphere to produce the average flexion over the 
whole sphere F. Taking advantage of the symmetry between the northern and southern hemisphere 
we can integrate only over the northern hemisphere (where dA = 2tt cos tpdcp) , yielding: 

F = (|/|) (13) 
/ 2tt cos 4> tan 4>^d(j) 
2~t~t 

2 

TT 

The flexion is less than that of the stereographic because of the 180° boundary cut along the 
longitude line at the international date line. A geodesic is a great circle on the globe, and if this is 
shown as a closed curve on the map that is always concave inward (the best possible case) it will 
always have a total rotation of the velocity vector of 360° and so will have an average integrated 
flexion of 1. If there is a boundary cut, the great circle does not have to close on the map (it has 
two loose ends at the boundary cut) and so need not completely rotate by 360°. 



3 Second Distortion: Skewness 

Acceleration in the direction parallel to the velocity vector of the truck a», causes the truck to 
increase its speed along the geodesic curve without causing any rotation. This causes skewness, 
because as the truck accelerates, it covers more distance on the map on one side of a point than on 
the other, so the point in question will not be at the center of the line segment of arc centered on 
that point on the sphere. 

Consider a segment of a meridian of longitude on the globe centered at 45° north latitude. 
Going from South to North along that geodesic in the Mercator map the truck is accelerating with 
ay > 0, because the scale factor is getting larger and larger the further north one goes, so as the 
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truck continues to cover equal arc length on the globe it covers larger and larger distances on the 

map. Thus, the center of the segment (at 45° latitude) is not centered on the segment on the map. 
We define the skewness: 

. = J (14) 
Taking the explicit case of the Mercator projection, we find: 

d y n k\ 

y y = # (15) 
i 



COS0 

and thus: 



tan</> 



COS(p 

so the skewness (for a vector pointed N-S) is simply: 

s = tan^ . (17) 

The skewness at 45° is 1, showing a lopsidedness toward the north. Given this relation, the 
skewness is positive (northward lopsidedness) in the northern hemisphere and negative (southward 
lopsidedness) in the southern hemisphere. 

Consider a geodesic through a point in the northern hemisphere tipped at an azimuth angle 
of 6 with respect to north. The only thing increasing the speed of the truck is the gradient of 
the scale factor as one moves northward, so the amplitude of the parallel acceleration is equal to 
the maximum acceleration (obtained going straight north) times cos 9. To get the average of the 
absolute value of the skewness for all geodesies through that point at all random angles 9, one 
simply integrates over 9: 

, . fn^ 2 cos 9d9 , , 

<s> = tan^P (18) 

7T/2 

,,2 

= I tan (j>\ — 

7T 

As with the flexion, we can integrate this over all points on the sphere to produce the average 
skewness over the whole sphere S. Similarly, we find: 

S=l (19) 

Notice that this is exactly the same value as the average flexion, F, for the Mercator. We will 
find that for conformal projections, the average absolute value of the skewness and flexion at a 
given point and over the whole globe are always equal. (This is only true for conformal projections, 
for general projections the skewness and flexion can be different, as illustrated by the gnomonic 
projection which has zero flexion but non-zero skewness.) 

In the Mercator projection, at the equator, the skewness is zero, as we would expect from 
symmetry considerations. Indeed, because any geodesic crossing equator has a symmetric shape 
in the northern and southern hemisphere, the skewness s = for any geodesic line evaluated at a 
point on the equator. Likewise, the flexion is zero for any geodesic line evaluated at a point on the 
equator. So the Mercator map has perfect local shapes along the equator, uniform scale along the 
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equator (Tissot ellipses all equal size circles) and zero flexion and skewness along geodesies in any 
direction from points on the equator. 

While much of the analysis in this work specifically addresses distortions in maps of the earth, 
these effects must also be taken into account in other maps as well. One of us (JRG) has recently 
produced a conformal "map of the universe" (Gott et al. 2005) based on the logarithm map of 
the complex plane. The horizontal coordinate is the celestial longitude in radians, yielding a 360° 
panorama from left to right. The vertical coordinate is: 



where d is the distance, and r ffi is the radius of the earth. The distance scale goes inversely as the 
distance, allowing us to plot everything from satellites in low earth orbit to stars and galaxies, to 
the cosmic microwave background on one map. The map is conformal, having perfect local shapes. 
However, it does have flexion and skewness. Circles of constant radius from the earth are bent 
into straight lines for example, and a rocket going out from the earth at constant speed would be 
slowing down on the map. 

In Figure O note that the continental United States also appears lopsided on the Mercator map. 
The geographical center of the continental United States (which is in Kansas) appears in the lower 
half of the continental United States on the Mercator map because the scale factor on the map 
gets larger and larger the further north one goes on the Map. Thus, the continental United States 
is lopsided toward the north in the Mercator map. Flexion or bending is manifest on the map as 
a bending of geodesies on the map, and skewness or lopsidedness is manifest on the map as the 
midpoint of a geodesic line segment on the map not being at the midpoint of that geodesic arc as 
shown on the map. 

4 Goldberg-Gott Indicatrices 

We began this discussion with the virtues of the Tissot indicatrices. Likewise, we have produced 
a simple indicatrix to show the flexion and skewness in a map as well as the isotropy and area 
properties indicated by the Tissot indicatrices. We will refer to them as "Goldberg-Gott" indica- 
trices. These are constructed as follows. At a specific point on the map draw a circle on the globe 
of radius 12°, and then plot it on the map. Inside this circle, plot the north-south, and east- west 
geodesies through the central point on the map. This leaves a © symbol on the map. If the map 
were perfect, this would be a perfect circle and the cross arms would be perfectly straight, with 
their intersection at the center of the circle. 

We have produced such a Goldberg-Gott indicatrix located at the geographic center of the 
continental United States for using a Mercator projection in Figure El 

We note that our technique is similar to that of Tobler (1964), who projected circles of 15 degrees 
in radius on to maps. However, it differs in that Tobler's circles do not have their centers marked 
nor perpendicular geodesies drawn from their centers like our indicatrices, and so do not convey the 
information on flexion and skewness. But they do show the finite shape distortions of the circles 
themselves which our indicatrices also cover. Tobler includes tables showing for particular locations 
on the sinusoidal projection, the maximum and minimum scale radius and maximum difference of 
radial directions for circles of various radii on the sphere (which address the Tissot issues of size 
and isotropy on finite circles but not flexion and skewness). Tobler also calculated errors in miles 
and angular errors for 300 random triangles within spherical circles of various radii in different map 
projections, and considered area errors for such random triangles within land areas in various world 
map projections. 




(20) 
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Using our indicatrices, one can see in Fig. [3] that the north-south geodesic is straight, but that 
the east-west geodesic is bent downward. This shows dramatically the flexion in this region of the 
Mercator map. One can even read off the average value of the flexion by hand. Take a protractor 
and measure the tangent to the east-west geodesic at the two ends of the cross bar. Measure the 
difference in the angle orientation of the two. That gives the integrated flexion along 24° of the 
globe. Divide that angle difference by 24° and you will have the average value of the flexion along 
that curve. 

The skewness is also visible in that the center of the cross is below the center of the circle, 
showing the lopsidedness to the north. In fact, one can observe the skewness in any direction from 
the center by seeing how far off center the center of the cross is with respect to the center of the 
circle in different directions. For comparison, we have in Figure [4] shown the continental United 
States in a oblique Mercator projection where the east- west geodesic through the geographic center 
of the continental United States is now the equator of the Mercator projection. 

The flexion and skewness along the equator of a Mercator map are indeed zero, so the arms of 
the cross are now straight, and the circle is now nearly a perfect circle centered on the center of 
the cross. This gives a "straight on" view of the continental United States, that more accurately 
portrays its appearance on the globe. 

One can place the Goldberg-Gott indicatrices every 60° in longitude and every 30° in latitude 
on the globe to show how the flexion and skewness vary over the map. In Figs. [5} [271 we provide 
G-G indicatrix maps for a number of well-known projections. 

In fact, the Goldberg-Gott almost indicatrices can just replace the Tissot indicatrices because 
the shape and size of the oval in the Goldberg-Gott indicatrix is approximately the size and shape 
of the Tissot ellipse. The Tissot ellipse shows the [magnified] shape and size of an infinitesimal 
circle on the globe, the oval in the Goldberg-Gott indicatrix © shows the shape and size of a finite 
circle (radius 12°) on the map itself at correct scale. Thus, if the map is equal area, the Goldberg- 
Gott indicatrices, will all have equal area on the map. If the map is conformal the Goldberg-Gott 
indicatrices will all be nearly perfectly circular. If there is a 2:1 anisotropy in the Tissot ellipses in 
a given region the Goldberg-Gott indicatrices ovals will have that same 2:1 axis ratio. 

5 A differential geometry approach 

Thus far, we have defined the general properties of skewness and flexion, given a few analytic results 
for particular map projections, and given a graphic approach for describing and interpreting flexion 
and skewness on maps. In this section, we approach the matter somewhat differently, and produce 
general analytic results for all projections as well as a prescription for measuring the flexion and 
skewness analytically. 

5.1 Coordinate Transforms 

Let's consider a spherical globe with coordinates: 



Note that here and throughout, we will use x a to refer to coordinates in the globe frame, and x' 
to refer to coordinates in the map frame. 
On the globe, the metric is: 




(21) 




(22) 
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such that, as always, the distance between two points can be expressed as: 



dl 2 = dx a dx b g-- b (23) 
Now, consider an arbitrary 2-d coordinate transformation: 

( x\<l>,\) \ 
~ { x 2 {^\) J 

Of course, from this definition, we may easily compute a local transformation matrix: 




(24) 



The inverse matrix is A° = dx a /dx a . The metric in the map frame is: 

g ab = A% Alg-- b (26) 
From this, we may then compute the Christoffel Symbols of general relativity: 

Tfec = ^ (9eb,c + gec,b - 9bc,e) (27) 

where standard convention tells us to sum over identical indices in the upper and lower positions, 
and where a comma indicates a partial derivative with respect to a coordinate. 

In practice, actually computing the Christoffel symbols for an arbitrary projection is not simple. 
To do this analytically requires that we have an analytic form for the map inversion. However, 
we make available a numerical code to compute the Christoffel symbols for all map projections 
discussed in this work on our projections website (see below). 

5.2 Analytic forms of Flexion and Skewness 

The whole point of computing the Christoffel symbols is that we want to address a very simple 
question: How are large structures distorted when projected onto a map? Clearly to an observer 
on the globe, a straight line is easy to generate. Point in a particular direction, and start driving 
(assuming your car can drive on the ocean) with the steering wheel set straight ahead. Drive for a 
fixed distance in units of angles or radian. Record all points along the way. 

Geometers, of course, know this route as a geodesic, and if we consider r to represent a physical 
distance on the surface of the earth, then the geodesic equation may be expressed as: 

= ~n c u b u c (28) 

where 

dr a 

u a = . (29) 
en- 



Equation (|28|) describes the bending of straight lines on a particular map projection, and thus, 
if all of the Christoffel symbols could vanish, we would clearly have a perfect Cartesian map. Not 
possible, of course. 

But what is the physical significance of the Christoffel symbols? Since the lower indices are 
symmetric by inspection, there are 6 unique symbols. What do they mean? 
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5.2.1 Analytic Flexion 

We define a vector oriented in a particular direction: 

8-W = ( 21 ) (30) 

In reality, this is not the unit vector, since: 

\u\ 2 = g ab u a u b £ 1 (31) 
Of course, we could define a true unit unit vector: 

u(9) = l(9)u(9) (32) 
where I is the "length" in grid coordinates of the unit vector. This has a value of: 

1(9) = 1 (33) 

J cos 2 6gu + sin 6g22 + 2 sin 9 cos 9gi 2 

Of course, it is clear that at any point on the map, the set of all u(9) represents an ellipse - the 
Tissot ellipse. 

We define the flexion in the following manner: Follow a particular geodesic a distance, dr 
(measured in, for example, radians). On the map, the geodesic will change direction by an angle 
d9. The ratio: 

f-Tr 

is the flexion. Note that for all polar projections, the equator will have a flexion of 1. 
In terms of the geodesic equation, the flexion can be expressed as: 

/(*) = W(* xu) 

= i(fl)(r^W-ry^ 2 ) 

(35) 

where u denotes a vector. 

Some light can be shed on equation (|35p if we expand the expression explicitly into trigonometric 
functions: 

f(0) = 1(9) [r 2 x cos 3 9+ (2T 2 12 - T^) cos 2 9 sinfl 

-(2T\ 2 - T 2 22 ) sin 2 9 cos 9 - T\ 2 sin 3 6>1 (36) 
This immediately reveals an important idea. If: 

r?i = r^ 2 = o 
rli = 2r 2 2 

T 2 2 = 2T\ 2 (37) 

then f9) in all directions will be zero. That is, all geodesies will be straight lines in this projection. 
As we will see, this is true for the Gnomonic projection. 

Inspection of equation (j36[) shows it is clearly anti-symmetric about exchanges of 9 and —9, 
since 9 is not the normal vector to the geodesic, but rather runs along it. 
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5.2.2 Analytic Skewness 

Essentially, a skewness means that if you walk (initially) north (for example) for 1000 miles, or 
walk south for 1000 miles, you will cover different amounts of map coordinate. 
The skewness along a geodesic can be defined similarly to the flexion: 



dr 

= 1(6) {v l bc u b u c u l + T 2 bc u b u c u 2 ) (38) 
As with the flexion, we may expand these out explicitly: 

s (0) = 1(0) [r} 1 cos 3 6» + (2T\ 2 + r 2 1 )cos 2 <9sin6> 

+(2T 2 12 + n 22 ) sin cos 2 9 + T 2 22 sin 3 6>] (39) 
A map with no skewness will have: 

rh = r 2 2 = o 
r 2 i = -2T\ 2 

V\ 2 = -2T 2 l2 (40) 



Unlike the flexion, we know of no projections with zero skewness everywhere. 

5.3 Projections with straightforward analytic results 
5.3.1 The Gnomonic Projection 

The Gnomonic Projection is particularly interesting. It has a coordinate transformation: 

(x \ __ I cot 4> cos A \ 
y J \ cot sin A J 

which is directly invertible to yield: 



(41) 



V A J \ tan" 1 (|) 
The coordinate transformation is thus: 



ros ! \/ ] (42) 



A|= y ^x 2 +y 2 (l+x 2 +y 2 ) y/x 2 +y 2 (l+x 2 +y 2 ) J ^43^ 

We thus have the map metric: 



y 

x 2 +y 2 x 2 +y 2 



g ab = ( ^jSr ) (44) 

(l+x 2 +y 2 ) 2 (l+x 2 +y 2 )' 2 
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We can compute the Christoffel symbols in the normal way. We find: 

2;r 



r 



i 



r?i = o 



11 1 + x 2 + y 2 

r 1 y 

12 1 + x 2 + y 2 

r 2 _ x 

1 + x z + y z 

rl 2 = o 



^ = (45) 
1 + x z + y z 

This clearly satisfies the requirements of equations f)3T[) . but not (|40p . and thus, the Gnomonic 
projection produces straight, but skewed geodesies. 

5.3.2 Stereographic 

The Stereographic projection is conformal, and thus, all of the Tissot ellipses are circles. Does 
this mean there is no skewness in the projection? No, as we've already seen. The Stereographic 
projection has the coordinate transformation: 



^ x j _ ^ tan(vr/4 + (/>/2)cosA 
which, again, is directly invertible to yield: 



y J y tan(7r/4 + (/>/2)sinA 



(46) 



<M - / f -2tan" 1 (^ 2 + y 2 ) 
X ~[ tau-i (f) 



(47) 



The coordinate transformation is: 



2x 2y 



V . 

x 2 +y 2 x 2 +y 2 



This can be used to compute the metric on the map: 



(l+x 2 +y 2 ) 2 ' 



ft*=| 4 ) (49) 



This clearly indicates that all Tissot ellipses will be circular. 

From these, of course, we can compute the Christoffel symbols: 

2x 



^22 — — ~ ^12 — T~, 2 (50) 

1 + x z + y A 



^li — — r 2 2 — — r 12 — — =— ■ 2 (^-*-) 

It is clear from inspection that geodesies are generally neither straight nor unskewed. 
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Moreover, it is clear that 1(9) is independent of orientation since the map is conformal. In 
general it can be shown to be: 

l = \(l + x 2 + y 2 ) . (52) 

In the stereographic projection, a circle of radius 12° on the globe is a perfect circle on the map 
but the center of the circle on the globe is not at the center of the circle on the globe (see Figure H]), 
and thus, there is skewness. 

6 Discussion 

6.1 Numerical Analysis of Standard Map Projections 

Not all projections produce such simple results. Thus, in general, we will want to compute the 
local flexion and skewness numerically. Our approach is as follows: For each projection we chose 
30,000 points selected randomly on the surface of a globe. For each of these points, we chose a 
random direction to start a geodesic. We follow that geodesic using small steps (dr ~ 10~ 5 rad) 
numerically, and use standard difference methods to compute the map velocity and acceleration 
along the geodesic. We are thus able to compute the metric and the Christoffel symbols (and thus 
the flexion and skewness) directly. We make our IDL (Interactive Data Language) code available 
to the interested reader at our projection webpage (see below). Likewise, we also do a distance 
test, in which pairs of points, (i,j), are chosen at random and the distance is measured both on the 
globe and on the map. 

This is a somewhat different perspective than simply inspecting the Goldberg-Gott indicatrices 
at a few locations, since we are now doing a uniform sample over the surface of the globe, rather 
than a uniform sampling over the map. When looking at the indicatrix map we can occasionally 
get a distorted view as to the quality of a particular projection. Some (like the Mercator) have 
relatively good fits over most of the globe, but the high latitudes can, in principle, be projected to 
infinite areas, and thus, the reader may erroneously think the Mercator infinitely bad. By sampling 
uniformly over the globe, we get a fair assessment of the overall quality of a particular projection. 

We define a number of fit parameters: I, corresponding to errors in the local isotropy (zero for 
conformal projections), A, corresponding to errors in the Area (zero for equal area projections), F, 
corresponding to flexion (defined in the discussion of flexion, above), S, corresponding to skewness 
(also defined above), D, corresponding to distance errors, and B, corresponding to the average 
number of map boundary cuts crossed by the shortest geodesic connecting a random pair of points. 

I = RMs(hx^\ (53) 

A = RMS (In aA- (In a t bi)) (54) 
F = (\fi\) (55) 
S = (\ Si \) (56) 

D = RMS (in (57) 

B = g (58) 

where Oj and bi are the major and minor axes of the local Tissot ellipses of random point, i, (Xi) 
indicates the mean of property, X, and Lb is the total length of the boundary cuts. 



13 



Logarithmic errors have been used before by Kavrayskiy (1958). Bugayevskiy and Snyder (1995) 
adopted schemes that weighted isotropy errors, (a/b — l) 2 , and area errors, (ab — l) 2 , according to 
one's interest in isotropy and area. Airy (1861) had just given these terms equal weight. 



Projection 


I 


A 


F 


S 


D 


B 


Non Interrupted Projections 














Azimuthal Equidistant 


0.87 


0.60 


1.0 


0.57 


0.356 





Gott-Mugnolo Azimuthal 


1.2 


0.20 


1.0 


0.59 


0.341 





Lambert Azimuthal 


1.4 





1.0 


2.1 


0.343 





Stereographic 





2.0 


1.0 


1.0 


0.714 





1 180 deg. Boundary Cut 














Breisemeister 


0.79 





0.81 


0.42 


0.372 


0.25 


Eckert IV 


0.70 





0.75 


0.55 


0.390 


0.25 


Eckert VI 


0.73 





0.82 


0.61 


0.385 


0.25 


Equirectangular 


0.51 


0.41 


0.64 


0.60 


0.449 


0.25 


Gall-Peters 


0.82 





0.76 


0.69 


0.390 


0.25 


Gall Stereographic 


0.28 


0.54 


0.67 


0.52 


0.420 


0.25 


Gott Elliptical 


0.86 





0.85 


0.44 


0.365 


0.25 


Gott-Mugnolo Elliptical 


0.90 





0.82 


0.43 


0.348 


0.25 


Hammer 


0.81 





0.82 


0.46 


0.388 


0.25 


Hammer- Wagner 


0.687 





0.789 


0.518 


0.377 


0.25 


Kavrayskiy VII 


0.45 


0.31 


0.69 


0.41 


0.405 


0.25 


Lagrange 





0.73 


0.53 


0.53 


0.432 


0.25 


Lambert Conic 





1.0 


0.67 


0.67 


0.460 


0.25 


Mercator 





0.84 


0.64 


0.64 


0.440 


0.25 


Miller 


0.25 


0.61 


0.62 


0.60 


0.439 


0.25 


Mollweide 


0.76 





0.81 


0.54 


0.390 


0.25 


Polyconic 


0.79 


0.49 


0.92 


0.44 


0.364 


0.25 


Sinusoidal 


0.94 





0.84 


0.68 


0.407 


0.25 


Winkel-Tripel 


0.49 


0.22 


0.74 


0.34 


0.374 


0.25 


Winkel-Tripel (Times) 


0.48 


0.24 


0.71 


0.373 


0.39 


0.25 


1 360 deg. Boundary Cut 














Lambert Azimuthal (2 hemisphere) 


0.36 





0.52 


0.11 


0.432 


0.5 


Stereographic (2 hemisphere) 





0.39 


0.37 


0.37 


0.692 


0.5 


Multiple Boundary Cut Projections 














Gnomonic Cube 


0.22 


0.37 


0.12 


0.87 


0.43 


0.686 



Table 1: The errors in isotropy, area, flexion, skewness, distances , and boundary cuts for some 
standard projections. 

In Table 16.11 we compare a these measures for a number of standard projections, which, for 
fairness of comparison we divide into a number of categories. The Gott-Elliptical, The Gott- 
Mugnolo Elliptical, and the Gott-Mugnolo Azimuthal have been discussed in Gott, Mugnolo & 
Colley (2007) and Gott, et al. (2007), where they have been applied to the earth, Mars, the moon, 
and the Cosmic Microwave all sky map. 

First, we show projections which represent the complete globe without interrupts. These pro- 
jections are azimuthal and the average flexion over these maps is 1. 

Second, we show the set of whole earth projections with one 180° interrupt. These include 



14 



rectangular and elliptical projections. Note that among all of the complete projections with or 
1 180° interrupt, there are a number of "winners" with regards to performance for flexion and 
skewness. The Lagrange has the smallest flexion. The Winkel-Tripel has the smallest skewness. 
For all conformal projections, the skewness is equal to the flexion. 

Of all of the whole earth projections, the most accurate for distance measure between points is 
the Gott-Mugnolo Azimuthal, followed very closely by the Lambert Azimuthal (Gott, Mugnolo & 
Colley 2007). 

In the third and fourth groups, we show 2-hemisphere and other multiple cut projections, 
respectively. 

In the final group, we have a projection with multiple interrupts, the gnomonic cube, which 
is defined piecemeal. This is a particularly interesting projection since the gnomonic is locally 
flexion-free, but it is clear that geodesies will not trace out straight lines in the gnomonic cube cube 
map because they bend when they cross an edge between faces. The gnomonic cube is presented 
as a cross, so 5 edges are included in the map proper. Geodesies bend when they cross an edge in 
this laid out cross configuration. 

The above comparisons do not depend on how important each of the criteria are (i.e. what 
weighting to give each measure). Laskowski (1997ab) suggested a means of ranking very disparate 
maps. Though his weighting scheme is not unique, as a simple illustration of how this can be 
done, we will minimize the sum of the squares of all 6 parameters, normalized to their values in 
the equirect angular projection: 

Mi)^)M^) 2+ © 2+ (0 2 +(0 2 

Following Laskowksi (1997ab) we set the normalization constants equal to the values of these errors 
in the Equirectangular projection [x = A, y = <p): TV, = 0.51, N a = 0.41, Nf = 0.64, N s = 0.60, 
N d = 0.449, N b = 0.25. 

The projections with the lowest values of E £ are: 

1. Winkel-Tripel (4.5629) 

2. Winkel-Tripel (Times Atlas) (4.5687) 

3. Kavrayskiy VII (4.8390) 

4. Gall Stereographic (5.7582) 

5. Hammer- Wagner (5.7847) 

6. Eckert IV (5.8519) 

This approach is certainly not unique. One may take issue with the weighting of the individual 
parameters, the domain over which they are applied (the whole earth, as opposed to continents 
only, for example), or even how the parameters are computed, we present it as a simple example of 
how our results may be combined with previous studies of map projections. It is interesting that 
the "best" map, as selected by this criterion is the one already used by the National Geographic 
for it's whole world projection. It is also interesting to note that the Winkel-Tripel has especially 
low skewness. 

Another nice property of this weighting scheme is that it gives a low score to pathological 
projections. For example, a series of n gores (made using the polyconic projection) arranged in a 
sunflower pattern would approximate the Azimuthal equidistant projection in distance errors as n 
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became large but would have arbitrarily low values of I, A, F, and S. But if a boundary cut term B is 
included this term would blow up and save us from picking the bad subdivided map as better than 
the more visually pleasing projections described above. Stretching individual pixels at the edges of 
the gores to fill the gaps between them would eliminate the boundary cuts. However, it would also 
cause the skewness to blow up, preventing us from giving a good score to a bad projection. 

Interested readers may visit www.physics.drexel.edu/~goldberg/projections/ to down- 
load a free IDL code to measure the flexion, area, and other measures discussed in this paper. 
We have not done all known projections, but have covered ones that have available mathematical 
formulas and we thought likely to do well. 

6.2 Conclusions 

We have developed two new measures of curvature distortions found in maps of the earth. The 
Skewness and the Flexion can be used to identify particularly warped sections of a map, or to 
identify features (such as the U.S. Canada border), which appear as a straight line on some map 
projections but do not, in fact, follow geodesies. We have developed a new graphical tool, called 
the Goldberg-Gott indicatrices, which can be used to show these distortions, along with those in 
area and shape, at many points in a map, and have produced indicatrix maps for a number of 
popular projections. Finally, we have used these measures to produce global distortion measures 
for different projections. We have found that the Winkel-Tripel produces low distortion on most 
measures, and, in particular, has the lowest skewness of all projections in our sample with a 180° 
boundary cut. 
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Figure 1: The Indicatrix map for a Stereographic projection. 




Figure 3: A Mercator projection cutout of continental the United States. We have have put a 
Goldberg-Gott indicatrix (a circle of radius 12° with N-S and E-W geodesies from the central 
point) at the geographic center of the continental U.S. 




Figure 4: A oblique Mercator projection cutout of the United States. We have have put a Goldberg- 
Gott indicatrix at the geographic center of the continental U.S. Notice that this projection has no 
flexion or skewness at the center. 
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Figure 5: The Indicatrix map for an azimuthal equidistant projection. 





Figure 7: The Indicatrix map for an Eckert IV projection. 




Figure 8: The Indicatrix map for an Eckert VI projection. 
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Figure 9: The Indicatrix map for a Gall-Peters projection. 




Figure 10: The Indicatrix map for an Equirectangular projection. 
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Figure 13: The Indicatrix map 



for a Gott-Mugnolo projection. 





Figure 15: The Indicatrix map for a Gott Equal-Area Elliptical projection. 




Figure 16: The Indicatrix map for a Hammer projection. 
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Figure 18: The Indicatrix map for a Kavrayskiy VII projection. 




Figure 20: The Indicatrix map for a Lambert Azimuthal projection. 
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Figure 24: The Indicatrix map for a Polyconic projection. 




Figure 25: The Indicatrix map for a Sinusoidal projection. 
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